Deconfined criticality in the frustrated Heisenberg honeycomb antiferromagnet 
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Using the density-matrix renormalization group, we determine tlie phase diagram of the spin 1/2 
Heisenberg antiferromagnet on a honeycomb lattice with a nearest neighbor interaction Ji and a 
frustrating, next-neighbor exchange J2. As frustration increases, the ground state exhibits Neel, 
plaquette and dimer orders, with critical points at J2/J1 = 0.22 and 0.35. We observe that both 
the spin gap and the corresponding order parameters vanish continuously at both the critical points, 
indicating the presence of deconfined quantum criticality. 



Introduction Models of frustrated magnetism on the 
honeycomb lattice have lately received tremendous inter- 
est. This interest stems from sign-problem- free Quantum 
Monte Carlo (QMC) studies which have established the 
presence of a spin liquid phase in the honeycomb Hub- 
bard model [H- Approaching from the strong coupling 
side, the physics at intermediate values of the Hubbard 
interaction U , for which the novel spin liquid phase has 
been found, can be described by the spin 1/2 Heisenberg 
model characterized by an antifcrromagnctic interaction 
Ji between neighboring spins and a frustrating, next- 
nearest neighbor exchange J2. When the frustration is 
small and J2 weak, the well-known Neel ordered state is 
stable, but at a critical value of a = J2/J1 it gives way 
to another, possibly liquid, phase. While all studies so 
far agree upon the presence of a phase transition, the 
nature of this intermediate phase that is reached by the 
transition out of the Neel state is heavily debated. The 
intermediate phase has been identified as a Z2 spin liq- 
uid by some [2-141 ^'iid as a plaquette-Resonating Valence 
Bond (pRVB) state, breaking translational symmetry, by 
others [5|-l7|. A recent variational calculation argues in- 
stead that the intermediate state does not have plaquette 
order Q . Upon further increasing the frustration param- 
eter a, a second transition takes place into a ground state 
that breaks lattice rotational symmetry but may or may 
not have magnetic order. 

We analyze this complex situation by formulating and 
answering four succinct fundamental questions on the 
Ji — J2 honeycomb Heisenberg model: (i) As to the Neel 
state: do quantum fluctuations tend to stabilize or de- 
stroy it? In other words, does Neel order vanish above 
or below the classical threshold of a = 1/6? (ii) What 
is the nature of the intermediate state? Is it a liquid 
state or docs it have plaquette order? (iii) What is the 
ground state for large a? Does it have magnetic order? 
(iv) What is the nature of the two phase transitions? Do 
the order parameters develop discontinuously or contin- 
uously across the quantum critical points? 

We use nominally-exact two-dimensional density- 
matrix renormalization group (DMRG) calculations to 
settle these issues and establish that: (i) Neel order is 
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FIG. 1. Phase diagram of the spin 1/2 Heisenberg antifer- 
romagnet on a honeycomb lattice with a nearest neighbor 
interaction Ji and a frustrating, next-neighbor exchange J2 
as obtained from DMRG. 



stabilized beyond the classical limit, up to a^i — 0.22 
(ii) the intermediate state has weak plaquette order with 
/-wave symmetry, and (iii) for ac2 > 0.35, the ground 
state has dimer order and breaks lattice rotational sym- 
metry. These results are summarized in the phase di- 
agram shown in Fig. [T] Moreover, we find that within 
numerical precision, (iv) both the spin gap and the rele- 
vant order parameters vanish continuously, at both crit- 
ical points Qfci and ac2- This implies that even if two 
different symmetries are broken on either side of ac, the 
transition is not first-order, as one would expect from 
a Ginzburg-Landau-type theory. Having two second- 
order transitions between the Neel, plaquette and dimer 
phases, implies that the critical theory for these tran- 
sitions is unusual and is not described in terms of the 
order parameter fields of either phase. It indicates in- 
stead the presence of two deconfined quantum critical 
points 

Frustrated honeycomb Heisenberg model The Hamil- 
tonian corresponding to the Ji — J2 Heisenberg model on 
a honeycomb lattice is 



H 



J, 



(1) 



where (ij) and {{ij)) denote nearest neighbor and next- 
neighbor sites i and j, respectively, and a ~ J2/J1 pa- 
rameterizes the strength of the frustration. We consider 
antiferromagnetic coupling: Ji , J2 and a arc all positive. 
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The model is well understood in the classical limit: at the 
critical value of a = 1/6, Neel order gives way to a spiral 
state with interesting order-by-disorder physics [ill Ell ■ 
However, in the extreme quantum limit of S = 1/2, 
the phase diagram is not well established We use 

DMRG to resolve this issue. 

Method Our DMRG is truly two-dimensional - we 
consider clusters with various geometries chosen to be 
conducive to various ordering patterns. It is well known 
that one can lift the degeneracy of wave functions by tak- 
ing some or all edges to be open. We use appropriate edge 
geometries as weak perturbing fields to induce symmetry 
breaking in the ground state. By performing measure- 
ments in the center of the cluster, one can estimate the 
order parameter induced by the edge geometry. Upon 
systematically increasing the size of the system, the ef- 
fect of the edges becomes progressively weaker and thus, 
by scaling to the thermodynamic limit, we can obtain the 
value of the order parameter in the ground state. In all 
cases, we have obtained smooth finite size scaling which 
indicates that our results exhibit a steady convergence to 
the thermodynamic limit. 

As described below, we have used a variety of cluster 
geometries appropriate for each phase. Note that the 
performance of DMRG calculation is equally stable for 
any ordered phase at a < 0(1). We study several cluster 
sizes with total number of sites up to 96 and keep up to 
6000 density-matrix eigenstates in the renormalization 
procedure. We perform ^ 10 sweeps until the ground- 
state energy converges within an error of ^ 10~^ Ji. All 
quantities calculated in this letter have been extrapolated 
to the limit n — > oo, where n is the number of retained 
eigenstates. 

Quantum stabilitization of Neel order We first de- 
termine the value of a at which Neel order vanishes 
and establish the role of quantum fluctuations in this 
process. Naively, one expects quantum fluctuations to 
destabilize Neel order for S = 1/2, thereby pushing the 
ad to a value below 1/6. On the other hand, as the 
Neel state is collinear, quantum fluctuations may pre- 
fer the Neel state over a competing spiral phase and 
push ad above 1/6. Even though various approaches 
have been used to resolve this issue, a consistent pic- 
ture has not emerged so far. Calculations which sup- 
port the hypothesis that ad < 1/6 include linear spin 
wave theory 13| , one-loop renormalization group study of 
the non- linear sigma model |14l |. functional renormaliza- 
tion group analysis 15|, and a Variational Monte Carlo 
(VMC) approach using RVB and Huse-Elser wavefunc- 
tions [J]. On the other hand, approaches which support 
the ad > 1/6 hypothesis include exact diagonalization 
(ED) 0, 0] , Schwinger boson mean- field theory [3 , se- 
ries expansions [l^], coupled-cluster calculations [3| and 
a VMC calculation using entangled plaquette states 

The DMRG results presented in Fig. [5] conclusively 
establish that quantum fluctuations stabilize Neel order 




FIG. 2. Finite size scaling of Neel order parameter. ( a) 
Diamond cluster with L — 3. (b) Hexagonal cluster with L = 
3. (c-d) Finite size scaling of Neel order paramater defined in 
Eq. [2] for diamond and hexagonal clusters, (e) Scaled Neel 
order parameter as a function of a = J2/J1 for diamond (red, 
closed circles) and hexagonal (blue, open circles) clusters. 



beyond the classical regime of stability. We have used two 
cluster geometries - diamond and hexagonal [Fig. [2ji,b] . 
One should be aware that periodic boundary conditions 
in some direction artificially enhance or diminish Neel 
correlations due to short range periodicity. This finite- 
size effect decays only slowly with increasing cluster size. 
To circumvent this issue, we keep all edges of the clusters 
open and measure the following order parameter as a 
function of a: 



711"^ [N] 



(2) 



As shown in Figj^l this quantity shows good finite-size 
scaling with terms proportional to 1/L and 1/L^, where 
L is the linear extent of the system. In the unfrustrated 
situation (a = 0), the staggered moment m in the ther- 
modynamic limit comes out to be 0.2857± 0.039 which is 
consistent with previously estimated values of 0.2677(6) 
and 0.270 obtained from QMC [l3| and ED respec- 
tively. As a increases, the obtained value of the Neel or- 
der parameter steadily decreases. At the critical value of 
ad ^ 0.22, we observe that Neel order vanishes in a con- 
tinuous transition as shown in FiglS^. Both diamond and 
hexagonal cluster geometries give the same value of ad , 
which signals the robustness of our result. Thus, quan- 
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FIG. 3. (a): Cluster geometry used to establish the presence 
of plaquette order. (b,c) Finite size scaling of spin gap and 
{Pcentrai) ~ a Hieasure of pRVB amplitude, (d) Spin gap and 
(Pcentrai) in the thermodynamic limit. 



turn fluctuations stabilize Nccl order significantly beyond 
the classical threshold. 

Non-linear spin wave analysis The excitations of the 
Neel state are well captured by spin wave theory, which 
treats quantum fluctuations using an expansion in pow- 
ers of 5". Linear spin wave theory with 0{S^) terms gives 
ad ^0.11 which is below the classical threshold. To 
reconcile this with the observed DMRG phase boundary, 
we take into account the quartic spin wave interaction 
terms of order 0{S^). Wc treat the interactions at mean- 
field level (for details, see Supplementary Material) and 
observe that the Hartrec Fock parameters merely renor- 
malize the strength of the Ji and J2 couplings. This 
effectively scales the frustration parameter a = J2/J1 
down so that the Neel state only becomes unstable be- 
yond a ^ 0.214. The quartic terms thereby provide a sig- 
nificant correction to the critical frustration ratio. The 
precise value of ad may depend upon further corrections 
beyond quartic order. Nevertheless, non-linear spin wave 
analysis confirms the strong tendency for quantum fluc- 
tuations to stabilize Neel order beyond the classical limit. 

Intermediate plaquette phase We observe the pres- 
ence of an intermediate plaquette-RVB (pRVB) phase, as 
suggested previously [iBlll, for 0.22 < a < 0.35. This 
state consists of a \/3 x arrangement of plaquettes as 
shown in Fig. [1] - each shaded plaquette is in an anti- 
symmetric combination of the two Kekule singlet covers. 
To test for plaquette order in the ground state, we choose 
the cluster geometry shown in Fig. [3^ which favors pla- 
quette order (this also favors columnar dimer order (l9|, 
but we have explicitly checked that it order does not oc- 
cur). This choice of boundary conditions acts as a weak 
field which induces plaquette ordering as shown by the 
shaded hexagons in Fig. [3^. To determine the pRVB or- 
der parameter, we first define the two single-plaquette 
states I a) and |6) - the two Kekule singlet covers of a sin- 



gle hexagon. The /-wave, antisymmetric, pRVB wave- 
function is given by |— ) ~ |a) — |6), upto a normalization 
constant. The order parameter corresponding to pRVB 
order is the projection onto the antisymmetric wavef unc- 
tion: OpRVB = |— )(— I acting on a shaded plaquette in 
Fig. [3^. We use the closely related plaquette- flip operator 
which flips the two Kekule covers: 



P 



\a){b\-\b){a\ 



(3) 



If the plaquette is in the pure |— ) state, this operator 
has expectation value 5 /4 (details in Supplementary Ma- 
terial). For the case of s-wave pRVB order, this expecta- 
tion value would be negative. 

To determine the strength of the pRVB ordering at the 
cluster center, we define (Pcentrai) as the average of (P) 
over three plaquette-ordering hexagons at the center of 
the system. As seen from our cluster geometry in Fig. |3^, 
one cannot always identify a single central plaquette for 
a given L. But we can always identify a central triad of 
plaquettes. Finite size scaling of (Pcentrai) provides the 
strength of pRVB order in the limit of infinite system 
size. Consistent with /-wave pRVB order, this expecta- 
tion value is positive for 0.22 < a < 0.35. Fig. [St shows 
the finite size scaling of (Pcentrai) which indeed scales to 
a positive value in thermodynamic limit. Also we find 
a finite spin gap that is consistent with \/3 x -v/S pla- 
quette ordering. We note, however, that strong quantum 
fluctuations reduce the amplitude of plaquette ordering: 
(Pcentrai) rcachcs a maximum value of ^ 0.43 compared 
to 5/4 for the case of pure pRVB order. The strength of 
pRVB order can also be characterized by p, the ampli- 
tude of the projection onto the |— ) plaquette wavefunc- 
tion as defined in Ref. [23|. For decoupled hexagons in 
the regime < a < 0.5, there is perfect pRVB order with 
p = 1. Our DMRG results indicate that in the honey- 
comb Ji — J2 model, pRVB order is strongly affected by 
quantum fluctuations and reduced to p < 0.43. To con- 
firm the existence of pRVB order, wc also measure the 
spin gap which is by definition the energy difference be- 
tween the first triplet excited state and the singlet ground 
state, 

A{L) = EiiL) ~ Eq{L), A= lim A(L), (4) 

L— ^00 

where En{L) is the n-th eingencnergy (n ~ corresponds 
to the ground state) of the system size L. The scaling 
analysis of the finite-size data is shown in Fig. and 
the results extrapolated to the thermodynamic limit are 
plotted in Fig. [3ji. We observe that the gap is finite only 
in the region of positive (Pcentrai)- 

Dimer phase At larger values of J2, the presence of 
a dimer state which breaks lattice rotational symmetry 
has been proposed previously [HI, 



181. This state has 



been variously called the staggered- Valence Bond Solid 
(s-VBS) or the Nematic VBS state in literature. We es- 
tablish that this state occurs in the phase diagram for 
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FIG. 4. Left: finite size scafing of (R), order parameter corre- 
sponding to lattice rotational symmetry breaking. For finite 



size scaling, we first take Lx ^ oo and then Ly 



We 



show data points obtained by Lx scaling for different fixed 
Ly values. The points connected by the line are final values 
obtained from Ly scaling. Inset: cluster geometry to detect 
dimer order. We enforce periodic boundary conditions along 
Ly and open boundary conditions along Ly. Right: The order 
parameter and spin gap (also obtained by similar finite size 
scaling) as a function of a = J2/J1. 



a > 0.35 using the cluster geometry in Fig|3^. We use 
open boundary conditions in the x direction and periodic 
boundary conditions along y, thus breaking the degener- 
acy associated with threefold lattice rotational symmetry. 
The cluster favors bond ordering with horizontal dimers 
as shown in FigH^. We first measure the breaking of lat- 
tice rotational symmetry by evaluating the expectation 
value of 



i? — Sa • Sr — Sr • Sr 



(5) 



where the sites A, B and C are chosen close to the center 
of the system (see Fig. \^). We determine the expecta- 
tion (R) while systematically increasing system size. If 
the true ground state breaks lattice rotational symme- 
try, we expect this quantity to scale to a non-zero value 
in the thermodynamic limit. For finite size scaling, we 
first take ^ 00 followed by Lj, — )■ 00. This sequence 
of limits ensures that there is no degeneracy arising from 
lattice rotations. We obtain smooth finite size scaling by 
restricting ourselves to even values of Ly, as shown in 
FigUl Including odd Ly values leads to small oscillations 
preventing smooth scaling. 

FigHb shows that (R) scales to a non-zero value for 
a > 0.35, clearly establishing broken lattice rotational 
symmetry in the ground state. However, this is consis- 
tent with two ground state candidates - dimer order or 
magnetic stripe order 0]. To distinguish between these 
two, we measure the spin gap. The finite size scaling for 
spin gap is shown in Fig. HJj. The error bars shown in 
Fig. |4)d are associated with the choice of to fit the data 
points. For 0.35 ^ a < 0.6, the spin gap scales to a 
non-zero value robustly. For a > 0.7, it is not possible 
to determine reliably whether the spin gap closes. The 
non-zero spin gap clearly indicates dimer order and rules 
out the presence of broken spin rotation invariance. 



Nature of phase transitions We have clearly demon- 
strated the presence of Neel, plaquette and dimer orders. 
Naively, one expects first-order quantum phase transi- 
tions (QPTs) between these phases as they break differ- 
ent symmetries. Our DMRG results, however, evidence 
a continuous transition out of the Neel phase: as can 
be seen from Fig. [2l the Neel order parameter vanishes 
continuously at ad = 0.22. This implies the presence 
of an exotic deconfined QPT Q. Approaching from the 
pRVB side, the quantum field theory governing this de- 
confined transition must involve spinous coupled to vor- 
tices in the pRVB order parameter [l^. This is an ex- 
citing proposition as a deconfined QPT in a model with 
realistic Heisenberg interactions has not been identified 
before. Surprisingly, DMRG results suggest that also the 
plaquette-dimer transition is continuous. As seen from 
Fig. [3ji and Fig. |4l at ac2 = 0.35, there is no evidence 
for either plaquette ordering or a breaking of lattice ro- 
tational symmetry. More detailed work will be needed 
to study the vicinity of these transitions, to extract crit- 
ical exponents and to rule out weak first-order behavior 
or the presence of a different small intervening phase. 
If there is indeed a continuous transition between dimer 
and plaquette phases, it would be yet another Landau- 
forbidden QPT within the same model. The field theory 
corresponding to this transition would be of immense in- 
terest. 

We thank I. Rousochatzakis for many useful discus- 
sions. 

Note added: during the preparation of this manuscript 
a DMRG study in Ref. [2l| reported a similar sequence of 
Neel, plaquette and dimer order as well as the continuous 
nature of the transition out of the Neel phase. 
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Non-linear spin wave theory 

To analyse the excitations about the Neel state, we 
follow the spin wave formalism of Ref. [ij. Using the 
Holstein Primakoff representation, the Hamiltoniann is 
expanded in powers of the spin length 5*. Ultimately 
however, we will set S = 1/2. The classical energy of the 
Neel state is given by terms proportional to 5^ 



ci 



3 J, 



where N is the total number of spins, 
correction, of order S, is given by 



^9" = Xl( b' 
















Ak J 





(1) 

The quantum 



(2) 



where 

Ak = S* 3 Ji — 6J2 + 2 J2{cos ka + cos kb + cos(fca + kb)} , 
Bk = -5Ji7k. 

We have defined 7k = e**' *, where (5's are the nearest 
neighbour vectors. The quantities ka and fcf, are compo- 
nents of momentum along two primitive lattice vectors of 
the triangular Bravais lattice. The operator aj^. (bl^) cre- 
ates a spin wave excitation on the A (B) sublattice. This 
Hamiltonian matrix can be diagonalized by a bosonic Bo- 
goliubov transformation with the eigenvalue 



(3) 



For J2 > Ji/6, the spin wave energy Ak becomes complex 
near the F point indicating that Neel order is unstable. 
We next include the quartic corrections arising from spin 
wave interactions by retaining terms of order 0{S^^^). 
There are no cubic terms. The interaction terms propor- 
tional to Ji are given by 



i,5 

+ajblb'lbi + alala^bl 



The index j stands for i + S. The quartic terms propor- 
tional to J2 are given by 

Hj2{0{S")) = -^X! [o'iO'lnO'l„a7n+ala^a^al„ 

b). (5) 



t t t t 



Aajaialj^ajn + (a 



The index m stands for i + rj, where 77 is a next-nearest 
neighbour vector. 

We treat these terms using the Hartree Fock approach. 
Using Wick's theorem, we replace bilinears with their 
expectation values and ignore the remaining normal or- 
dered quartic interaction piece. We take only the follow- 
ing bilinears to have non-zero expectation values: 



{a} at) 
(aibi+s) 



P, 

: h. 



= n, 



(6) 
(7) 
(8) 



These are the only bilinears which have non-zero expec- 
tation values within the quadratic theory. This choice 
of order parameters leads to a self-consistent theory that 
does not induce extra bilinears with non-zero values. The 
quantities n and h, being expectation values of Hermi- 
tian operators, are real. We take p to be real, as it is real 
within the quadratic theory. Using the symmetries of the 
underlying Neel state, we take these three quantitites to 
be independent of position and choice of neighbour (6, rj) . 
The interaction terms can be decoupled as 

Hj,{0{S")) = JiY,[Hp-n}{alau + blb^) 

k 

+ {n-p}7k(a-k&k + ak^Lk) • 
Hj,{0{S")) = J2^ [6{n~h}{alau) 

k 

- {n- /i}/Xk(aj^ak) +{a^b). (9) 

Here, /ik = J^r^^^Pi^^ ' v)- With this decoupling, these 
terms enter the quadratic Hamiltonian in Eql21 We have 



A + 3 Ji{p - n} - 6J2{h - n} + J2{h - nj^k 



and 



B ^ B + Ji{n ~ p}jk- 



Clearly, the Hartree Fock decouplings merely renor- 
malize the exchange couplings Ji and J2. We have 



Ml + {p-n}/S); 



J2 



J2{l + {h-n}/S). 



We obtain the Hartree Fock parameters n, h and p self- 
consistently. For every 'bare' value of J2/J1, we obtain 
a renormalized value of J2/J1 as plotted in Fig. [T] Figl2] 
plots the obtained value of n. When n ^ S, the Neel 
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FIG. 1. EfTective J2/J1 obtained after Hartree-Fock treat- 
ment of interactions. Note that the effective J2/J1 lies only 
approaches the instability threshold 1/6 when the bare ratio 
~ 0.95. 



0.2 0.3 



FIG. 2. Plot of the order parameter n. When n ~ 5, 
Neel is expected to become unstable to quantum fluctua- 
tions. For S = 1/2 (shown as a dotted line), this happens 

at (J2/Jl)i,are ~ 2.14. 



lattice 0. Thus, our value of critical J2/J1 is close to 
Schwinger Boson mean- field result of Ref. Q . 



Plaquette operators 

On a single plaquette, we denote the two Kekule sin- 
glet covers as \a) and \b). We take these states to 
be normalized. They are however not orthogonal with 
{a\b) = —1/4, upto a phase that can be gauged away. We 
denote symmetric (s-wave) and antisymmetric (f-wave) 
combinations of these covers as 



(l«> + l^'»,l+> = 



5(l«>-|^'». 



(10) 



With this definition, we have (+|+) = (^|~) ~ 1 and 
(+|— ) = 0. On an isolated plaquette, this operator takes 
the expectation values: 

{a\P\a) = 1/2, {h\P\b) =1/2 (11) 

(+|P|+) = -3/4, (12) 

(+l^'|-)-0, (13) 

(-|P|-)=5/4. (14) 

Our cluster supports pRVB order. With pure pRVB 
ordering, the central plaquette must be in the |— ) state. 
Wc find a positive expectation value for Pcentrai which 
supports the hypothesis that the central plaquette is in 
the |— ). However, the expectation value is less than 
5/4 which indicates that quantum fluctuations reduce 
the strength of pRVB order. Our cluster geometry can 
also accommodate c-VBS order, in which case the central 
plaquette would be in a pure \a) or \b) state. We have 
explicitly checked that this does not occur. 



moment is rcnormalizcd to zero and Neel order is ex- 
pected to become unstable to quantum fluctuations. For 
S = 1/2, this happens for {J2/ Ji)bare ~ 0.214. 

Our self-consistency equations are equivalent to a par- 
ticular formulation of Schwinger Boson mean-field the- 
ory. This connection between two very different methods 
has been pointed out earlier for the case of the square 
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